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In this paper, we numerically study the forward correlations between the lepton-pair and asso- 
ciated hadrons in Drell-Yan process in pA collisions. Using the present knowledge of the dipole 
gluon distribution from the modified Golec-Biernat-Wiisthoff model and from the solution of the 
Balitsky-Kovchegov evolution equation, we are able to compute and predict the forward correlations 
between the lepton-pair and associated hadron in Drell-Yan process at RHIC and LHC. Similar to 
the forward dihadron correlation in d-Au collisions measured at RHIC, the Drell-Yan type corre- 
lation also implies a strong suppression of the away side hadron at forward rapidity due to the 
strong interaction between the forward quark from the projectile proton and the gluon density from 
the target nucleus. Another feature of this process is that the correlation contains a double-peak 
structure in the away side, which makes it a unique observable. 



I. INTRODUCTION 

Over many years of both theoretical and experimental studies, Drell-Yan lcpton pair production has been considered 
as one of the most interesting processes in high energy physics. Especially, single inclusive Drell-Yan lepton pair 
production in pA collisions is a cleaner probe than high-pj^ hadrons since it has neither final state interactions nor 
fragmentation effects. Therefore, it provides a unique opportunity to study the parton distributions, especially the 
unintegrated gluon distributions [TH3] , in hadrons. Furthermore, according to the high energy factorization, the 
hard photon-hadron and Drell-Yan lepton pair-hadron (virtual photon-hadron) correlations in p^4 collisions serve an 
important role in probing the low-x dipole gluon distributions. 

In the small- a; regime [4], there exist two distinct unintegrated gluon distributions, namely the Weizsacker- Williams 
(WW) gluon distribution [SJ [5] and the dipole gluon distribution [7] . According to its operator definition, the WW 
gluon distribution corresponds to the conventional gluon distribution, which measures the gluon density in the light- 
cone gauge. Although the dipole gluon distribution has no such interpretation, it is a fundamental quantity which 
appears in many processes involving hadrons. The dijet (or dihadron) correlations in DIS [5l41Uj can help us directly 
measure the WW gluon distribution which has never been explored experimentally before, while the Drell-Yan lepton- 
pair-hadron correlation can directly probe the dipole gluon distribution at small-x. Since the virtual photon does not 
interact with the gluons inside target hadrons, the final state effects are absent. As a result, this correlation can be 
calculated exactly in the leading order for all angles between the virtual photon and associated jet or hadron [TO]. In 
contrast, due to the final state interactions, the so-called correlation limit or back-to-back limit has to be taken in 
order to arrive at a factorized formula for the dijet or dihadron processes. 

Previous studies |llH17j of the Drell-Yan processes in pp or pA collisions have focused primarily on the inclusive 
cross sections, but not on the correlations between the outgoing hadron and the Drell-Yan lepton pair 1 . As shown 
in Ref. [8141 0], the correlations can reveal more dynamical information about the dense nuclear target and yield the 
direct measurement of the dipole gluon distribution. 

The objective of this paper is to numerically study the correlations between the lepton-pair and associated hadron 
in the Drell-Yan process in pA collisions. This provides predictions for the correlations which can be measured at 
both RHIC and LHC Ql]. 

To parametrize the dipole gluon distribution, we have employed two approaches. In the first one we used the Golec- 
Biernat-Wiisthoff (GBW) model j2Q] together with the geometrical scaling [21] which can successfully describe all the 
low-cc DIS data. In deep inelastic scattering, the partonic cross section for a virtual photon scattering off a proton can 
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1 We also notice that there is a recent studv[T8] on the real prompt photon and hadron correlations in pA collisions. If we take the 
virtuality of the photon M 2 to be 0, and subtract the collinear divergence with proper scheme, we can also get the results for the real 
prompt photon and hadron correlations. 
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FIG. 1: The Drell-Yan scattering process, with several of the momentum variables used in our calculation labeled. p p and 
Pa are respectively the momenta of the proton and the nucleus from which the gluon was emitted. p 7 is the momentum of 
the virtual photon as reconstructed from the lepton pair, and p-n is the measured momentum of the pion. x p and x g are the 
longitudinal momentum fractions of the quark relative to the proton and the gluon relative to the nucleon, z = p^ /q + is the 
fraction of total momentum taken by the photon, and Z2 = pt /kg is the fraction of the momentum of the quark jet that is 
carried by the pion. 



be written in the color dipole model as a convolution of the photon wave function and the dipole-proton scattering 

(21 

amplitude. In the GBW model, the dipole scattering amplitude is parametrized as Sq-^ w (t^) = exp 



_ 

r_L being the transverse size of the qq dipole. The geometrical scaling property means that if one writes Q 2 {x) 

QsQ 



with 



2 (x/xq) x with Q s0 = lGeV, x = 3.04 x 1CT 3 and A = 0.288 [22], the DIS total cross section only depends 



on one single variable r = , instead of the variables x and Q 2 separately. Although the geometric scaling was 
originally found in DIS on proton targets, it can be generalized to pA collisions as well [23]. One can write the 
saturation momentum scale Q s a for target nuclei with mass number A as Q 2 sA (x) — Q 2 Q As c(b)(x / xq)~ x , where c(b) 
is the profile function which depends on the impact parameter b. This profile function is related to the centrality of 
the pA (or dA) collisions. Central collisions give large values of c(b), while the peripheral collisions correspond to 
small values of the profile function. 

In the second approach, we numerically solve the small-x evolution equation for the dipole amplitude, namely the 
Balitsky-Kovchegov (BK) equation [24] [25], using both fixed and running couplings. From these numerical solutions 
we then obtain the dipole gluon distribution in each case. We will discuss the implementation of this approach later 
in more detail. 

Comparing to the dihadron (or dijet) correlations in pA processes, the Drell-Yan type correlation contains a unique 
feature on the away side. This correlation has a double-peak structure in the away side while other correlation functions 
exhibit a single peak structure. This difference comes from the fact that the dipole gluon distribution and the cross 
section both vanish when the transverse momentum of the gluon goes to zero. This leads to a vanishing partonic cross 
section when the produced virtual photon and quark are completely back-to-back. With the fragmentation effect 
which turns the quark into a pion (tt°), this effect gets smeared a little bit and becomes less visible. Nevertheless, 
taking into account the fragmentation, we find that the correlation has a minimum at Aip = n and possesses a 
two-peak structure. 

The rest of the paper is organized as follows. In Sec. [Tl] we discuss the cross section of the Drell-Yan process with 
associated hadron in pA collisions and define the correlation function by dividing the single inclusive cross section 
of the DY process. We discuss the parametrization of the dipole gluon distributions in terms of different models in 
Sec. Ill The numerical results and further discussions on the future experiments are given in Sec. IV In Sec. [V] we 
conclude and summarize. 



II. DRELL-YAN LEPTON PAIR - HADRON CORRELATIONS 



In the small- a; regime, the Drell-Yan lepton pair production processes is dominated by the qg — >• qj* channel at the 
partonic level. Here the gluon g, in fact, can be regarded as a group of soft gluons from the dense nucleus target. The 
virtual photon then decays into a lepton pair and the quark fragments into a cluster of hadrons. We are interested 
in the correlation between the lepton pair and the hadron, say 7r°, which is a product of the fragmentation from the 
quark in the process pA — > £Itt°X. It is shown schematically in figure [l] 

The cross section AaP dv g > which measures the probability of producing a virtual photon with an invariant 
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mass M and a quark with large p± in the final state, has been calculated in [TU] for the pA —> "f*qX channel. After 
summing over different photon polarizations, the total production cross section can be cast into 2 



da pA ^ 1 * qX 
dy 7 dy,d 2 p 7 ±d 2 k 9± d 2 6 



2ir 2 



(l-z)F x (q ± ) 



[P 



2 2 



e 2 M ] [(P ± + zq ± 
-z 2 (l-z)M 2 



:: A/ 



(i) 



where e| f = (1 — z)M 2 , q± = p~ ( ± + k g j_ and Pj^ = (1 — z)p 7 j_ — z\i q ±. Here the normalized gluon distribution, 
Fx g (g±) — J f 2 l)2 e~ zq - L " r± S X 2 J (r±), is defined through the Fourier transform of the dipole amplitude. It is straight- 
forward to write down the corresponding cross section for the lepton pair production as follows 

d(J pA^i+r q x d(J pA^ q x 



dY 1 dY q d 2 p 1± d 2 k q± d 2 bd 2 M 2 3ttM 2 dK Y dY,d 2 p 7± d 2 k, ± d 2 & 



(2) 



Therefore, one just needs to include a factor of 9 ^ L< j^p- to the photon cross sections to get the lepton pair cross 
sections. The factor ^^jji eventually cancels out in the correlation function (J5j|. 

To express the cross section in terms of the kinematic variables of the pion (tt°) and leptons, which are actually 
detected in the final state, we need to include in the above formula the fragmentation function of the quark into a 
pion, D^o /f{%2, fJ-) to obtain 



dr 7 dy 7r d 2 p 7 _ L d 2 p 7r± d 2 6 



dz 



y 2 X! 



ftgJJ] G £ 

D 7 ,o /f (z 27 n)x p qf(x p , 



[1 + (1--) 2 



[P. 



2 2 



1 [(Pj j- ■<-< i 2 



AI 



2/1 _^n,f2 



C MJ 



z 2 {l-z)M 2 



Pi 



-M 



-i 2 - 



(P± +zq ± ) 2 + e 



71/ 



(3) 



where zhi = P 7 /pp~ and Zh2 = pt /Pp are the longitudinal momentum fraction of the virtual photon and final state 
7T° hadron, respectively. fi 2 is the factorization scale which is set to be Q 2 ^- The lower limit on the z 2 integral comes 
from requiring that x p < 1. 

The Drell-Yan lepton pair inclusive differential cross section can be obtained from the previous expression by 
integrating over the phase space of the final state quark, which makes it a function of the kinematic variables of the 



virtual photon only. By noting that dY 7 dY g 



dzdx p 
z(l-z)a 



-, one can write the inclusive Drell-Yan cross section as 



d(J pA^'X 

dY 7 d 2 p 7 j_d 2 6 



dz 

z 



d 2 q_L *S2 x p q f (x p , [i) - 



f 



2ir 2 



-F Xg (q±) 



[1 + (1--) 2 ] 



2 2 

z q± 



[p 2 ± + 4i] [(pi± - zc i±-) 2 + 4A 

i 



z z (l - z)M l 



P 2 ± 



1 AI 



(PjX ~ zqj_) 



-M 



(4) 



2 The rapidity is defined as respect to the center of mass frame of the scattering, which coincides with the lab frame when the energy of 
the proton projectile is the same as the energy per nucleon in the nucleus target. 
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This result is in agreement with the results in Ref. [T2lU5| . 

Integrating Q and Q over the magnitude of the transverse momenta with certain cutoff gives the cross sections 

as a function of rapidity and impact parameter only, dY d Y d 2 b and dY d 2 b • ^ se ^ ^ ne l° wer limits on these 
momentum integrals, we impose a transverse momentum cutoff p_L C ut which can be set to match the experimental 
acceptance. The upper limits are again determined by requiring that the longitudinal momentum fraction of the initial 
quark x p be less than 1. This would give the lower bound of zi integrals and upper bounds of momentum integrals. 
For the correlation measurement, either the virtual photon or the detected hadron can be the leading particle with 
larger transverse mass, while the other becomes the associated particle. In the code, we break the total contributions 

into two cases: (1) \Jp 2 ±_ + M 2 > Pq±', (2)^/p 7 j_ + ^ 2 < PqJ- which correspond to the region z > | and z < |, 
respectively. 

Finally, the ratio of these two cross sections defines the Drell-Yan pair-hadron correlation function C DY (A</>), 



2tt J - J dp 7 _ L p 7 _ L dp 7r _ L p 7r _ L 



t j .pA-y 7 * 7r° x 

dy 7 dy T d 2 p 7 j_d 2 p ff _Ld 2 6 



C DY (A</)) = . . , ( 5 ) 

d 2 p 7 _L- 



da pA^'X 

-dY 7 d 2 p 7± d 2 6 

Pj±>P±cut 

where A</> is defined as the azimuthal angle difference between the virtual photon and measured ir° . The correlation 
function expresses the probability density for a pion to be emitted at a specific azimuthal angle relative to the virtual 
photon (normally the rapidities of these two particles are chosen to be the same). When A0 ~ n, that is, when 
the virtual photon and ir° are almost back-to-back, the dominant contribution to the correlation function C DY (A0) 
comes from the low q± region of the gluon distribution F Xg (q±). In contrast, the behavior of the correlation function 
C DY (A0) at Acf> ~ or 2ir is mainly determined by the large transverse momentum part of the gluon distribution 
function F Xg (q±). 

In addition, as mentioned earlier, the partonic cross section vanishes at qj^ = as in Eq. (J3j) . The physical inter- 
pretation is that the emission of the virtual photon from the quark requires that the quark gets a kick or an impulse 
in the transverse direction. Since the small- a; gluons carry very little longitudinal momentum and the virtual photon 
does not interact with the dense small- a; gluons from the target nucleus, the quark must change its transverse mo- 
mentum during the multiple scattering with the gluons in order to radiate the photon. In other words, if the quark 
receives no impulse which change its momentum, the virtual photon will remain as part of the quark higher fock state 
wavefunction, and thus can not be produced in the final state. Due to this fact, we expect that A</> = n should be 
a local minimum of the correlation function. Furthermore, there are two symmetrical maxima lying near A</> = ir as 
q±F x (q±) increases with qj_ and reaches its maximum at q± — Q s . As a result, the correlation should have a unique 
double-peak structure on the away side. Note that the correlation function at A</> = n is no longer zero due to the 
pion fragmentation function. 

In terms of the numerical evaluation of the cross sections ^ and Q, we use the MSTW 2008 NLO parton 
distributions |26j for the parton distributions qf and the DSS fragmentation functions [27 for the fragmentation 
functions D^o / g . In the case of the d-Au collisions at RHIC, since the parton distributions are only given for protons 
in Ref. [26], we obtain the parton distribution for neutrons by using the isospin symmetry, which gives q u / p = qd/n 
and q d/p = q u/n . 

III. THE DIPOLE GLUON DISTRIBUTION 

The Drell-Yan cross section involves the so-called normalized dipole gluon distribution F Xg (q±), which characterizes 
the dense gluon distribution in this process in the nucleus in the small-x regime. As discussed in Ref. [5], this function 
is related to the well-known dipole gluon distribution xG^ as 

xG^(x,q ± ) = ^qlF Xg ( qi _), (6) 

where S± is the transverse area of the target nucleus, N c — 3 is the number of colors, and a s is the strong coupling 
constant. 

As mentioned previously we have chosen two approaches to parametrize the normalized dipole gluon distribution 
F x (<L_)j namely, the solution of the BK evolution equation [52] [5S] and the GBW model [20... F x is a Fourier 
transform of the color dipole amplitude in the coordinate space. The implementation of the solution to the leading 



5 



order BK equation for the phenomenological use is a little bit troublesome since it gives too fast energy evolution as 
compared to the experimental data from DIS. The GBW model describes the DIS data very well and is straightforward 
to use. However the GBW model fails to describe the proper behavior gluon distribution at large transverse momentum 
since it has an exponential drop off, while the pcrturbativc QCD gives a result which drops in terms of powers of q± . 
In practice, we use these two approaches for our numerical evaluation since they are complementary to each other. 
Our procedure is to fine-tune the behavior of the strong coupling constant a s or the running coupling a s (q) in the 
numerical solution of the BK equation until the solution has the same energy dependence and similar low q± behavior 
as what the GBW model gives. Then we believe that the solution of the BK equation should be able to produce the 
correct large q± part of the distribution and thus the correct angle correlation at A<p = 0, 2n. 



A. The BK Equation 



The Balitsky-Kovchegov (BK) equation governs the rapidity dependence of the color dipole scattering amplitude 
jV(r, b, Y) [211 [25]. In position space, we can write the equation in terms of S(r, b, Y) = 1 — Af(r, b, Y), as 



as(x- y ,b,y) 

dY 



(x-y) 



2tt (x-z) 2 (z-y) 



S(x-y,b,Y) 

-5| x z,b 



z - y 



,Y )S z-y,b- 



■ Y 



(7) 



The above form of the equation is valid in the leading logarithmic order in In 1/x. One can also write the BK equation 
in terms of the momentum-space dipole scattering amplitude, 



d 2 r 



. r N(r,Y) 



As shown in Rcf. 



in terms of the function <fi(k, Y), the BK equation reads 
d<t>(k,Y) 



dY 



= a s K (g> </>(k) - ot s (j) Z (k) 



where K represents the action of the BFKL kernel, which in the leading logarithmic approximation reads 

dfc' 2 rfc' 2 0(fc') - fcvo) k 2 4>{k) 



K (g> (j>{k) = 



V4fc' 4 + fc 4 



The Laplacian of <fr is the function F Xg that appears in the cross sections of section |TTJ 

i^ a (k,b) = i-V 2 0(k,b,r(x 9 )) +,5 2 (k) . 



(8) 



(9) 



(10) 



(11) 



Since the partonic cross section vanishes at k = 0, we can safely ignore the delta function in the above equation. 
Beyond this section, the b dependence is left implicit. 

The normalized dipole gluon distribution F x at small x extracted from the solution to the leading order BK 
equation yields energy dependence which is too fast as compared with phenomenology. In addition, the numerical 
solution to the full NLO BK equation [3S] is not yet ready for use in phenomenological study. In practice, we have 
chosen to modify the solution from the LO BK equation to obtain a realistic energy dependence by varying the 
behavior of the coupling constant. For a fixed coupling, we vary the value of the coupling constant directly; for the 
running coupling, we vary the value of Aqcd which controls the running behavior. 



B. The Golec-Biernat Wusthoff Model 



We can also parametrize the dipole gluon distributions by using the GBW model adapted for nuclear targets, which 
has both the nuclear enhancement and geometrical scaling [211 ED] feature. The geometrical scaling was found to be 
related to the traveling wave solutions [5TH55] of the Balitsky-Kovchegov evolution. 

For a dense proton target, the GBW model reads 

J d 2 bAA(r,b,y)=AA (l-e-^/ 4 ) , Q 2 (F) = Q^xfc™ = Q% (^j \ (12) 



() 



where Q 2 — lGeV 2 , A = 0.288 and x$ = 3.04 x 10 4 . As for a dense nuclear target with A nucleons inside, the 
saturation scale is subject to the nuclear enhancement which modifies the saturation momentum and gives Q 2 sA {x) — 
Q 2 a A^ c(b)(x I 'x )~ x . Here the profile function c(b) gives the centrality of the scattering. Using the GBW model, we 
can compute a simple analytic expression for F x . According to equations (|S[> and (11) above, we get 



4>(k\Y)= l -T[u. 

F Xg (k 2 ,Y) 



k 



o-k 2 /Qi A (Y) 



(13) 
(14) 



where r[0, a;] is the incomplete gamma function. 



C. The numerical solution of the BK equation 



We solve the BK equation numerically [25] with initial condition for <fi(k 2 ,Y) as in Eq. (13 1 at Yj n it = — lnx; n it, 
where we use Xi n it = 0.01. This value is chosen to coincide with the upper applicable limit of the GBW model. 
With a solution for <j), using F Xg (k 2 ,Y) — 9 a{\nky ^ ' wc tncn com P u te F Xg in the region Y > 5^ n it = 4.6 using 
a proper numerical differentiation scheme. The numerical solution to the BK equation for the gluon distribution is 
shown in figure [2] with fixed coupling, and in figure [4] for running coupling, as a function of the transverse momentum 
for different rapidities. For comparison also the analytic expression from the GBW model (14) is shown in figure [2] 
The BK solution possesses a larger tail in momentum as compared with the GBW model, which is due to the correct 
perturbative QCD limit. 



1. Fixed Coupling 

Since the solution to the leading order BK equation gives energy dependence which is too fast as compared to 
data (for inclusive DIS) if we use the normal value for a s , we treat the strong coupling constant a s as an adjustable 
parameter in the numerical computation. We need to find out the value of a s which corresponds to the kinematic 
regime probed by these collisions in order to make meaningful predictions. Theoretically, the value of A is predicted 
by various analytical approaches to the solution to the BK equation [3TI 153"] . The resulting saturation scale can be 
cast into 

Qt(Y) = ^^eJ-^±Y), (15) 



Y 2(1-70) V 1 ~ 7o 

where a s = ° , Ar ' , x(l) is the usual characteristic function the value of 70 is determined to be around 0.37. This 
function asymptotically approaches an exponential form, 

Q 2 (Y)^CeM~m- whereA=^^ = i^a s . (16) 

l-70 7T 

In terms of the numerical computation, as shown in figure [3lwe extract the saturation scale from the output of the 
BK evolution and employ the procedure as described in Ref . |28| . determining the value /c max of k which maximizes 
k(j>(k 2 ,Y) at each Y, and set the saturation scale to be proportional to fc max . It is straightforward to show that 
the assumption fc max oc Q s is valid for the GBW model. By matching the asymptotic behavior of the saturation 
momentum at large Y with the saturation momentum in the GBW model, we find that a s ~ 0.062 which is unusually 
small value for the strong coupling constant. Such small value of the strong coupling constant is an artifact of the large 
higher order corrections to the BK equation. If the higher order corrections and possibly resummation is included, 
the energy dependence extracted from the numerical solution of the BK equation should be close to the measured 
value with a reasonable value of a s (see, e.g., in Ref. [M]). 

With the fixed coupling, it seems that the BK evolution produces a transient spike in F x at low k 2 , which is 
primarily due to our choice of the initial condition which differs substantially from the asymptotic solution at low k 2 . 
Given the fact that our code samples F Xg within some extended region of q\ for a certain value of A<j>, this has very 
limited effect on our results. 
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FIG. 2: The gluon distribution from the GBW model (14 1, on top, and the output of the numerical BK evolution with fixed 



coupling (111, on bottom, for central p-Pb collisions. Each curve shows F Xg (k) for x g = exp(Yi n i t — n8Y), where n is the 



number that labels the curve and 8Y = ^ ln(10 8 ) w 0.04605 is the step size in rapidity used by the numerical integrator. The 
curve labeled is the initial condition. Note the transient "spike" at the low k 2 region of the solution to the BK evolution, and 
the significant enhancement at large k 2 relative to the GBW model. 



2. Running Coupling 

Incorporating running coupling corrections into the leading order BK evolution is known to give a reasonably 
accurate description of DIS data [3SH3B], and also to slow down the energy dependence of the saturation scale, which 
makes it a viable and perhaps preferable alternative to artificially tuning the value of a s . In this case, the saturation 
scale dependence has also been calculated in [31] to be 



«i -70) 1 VVW-ToJxho) 






a s = 0.2 




AgjcD = 0.0588 




GBW 




A^ CD = 0.0010 




a a = 0.062 



Y 

FIG. 3: This plot shows the peak of the momentum distribution, fe max = max k(f>(k 2 , Y), computed from the analytic formula 
for the GBW model, from the fixed coupling BK evolution for selected values of a s , and from the running coupling BK evolution 
with selected values of Aq Cd . For the BK evolution curves, the slope in the upper range of rapidities decreases as a s or Aq CD 
decreases, and the closest match to the slope of the GBW model curve at Y > 15 is achieved with a s — 0.062 for fixed coupling 
or Aq CD = 0.001 for running coupling. The jagged "steps" in the curves reflect the finite spacing of the momentum grid used 
in the evolution. 

with a running coupling of the form 



a s (k 2 ) 



f3lnk 2 /A 2 



(18) 



QCD 



where 



ll-2N f /N c 



l , £i = —2.338 ... is the rightmost zero of the Airy function, Nj = 3 is the number of quark 
flavors that are expected to contribute, and N c = 3 is the number of colors as before. Ref. [33] also includes a similar 
result. 

Aqcd- We use 



However, this form is not practical for the numerical implementation due to a Landau pole at k 
an implementation which shifts the squared momentum, 



a s (k 2 



1 



/Sin 



A 2 

QCD 



(19) 



with /Lt = 1 GcV. The shift term jj, 2 is a purely phenomenologically motivated change, which does not appear to 
significantly alter the behavior of Q S (Y) discussed below. 

In equation (191, we are treating Aq CD as an adjustable parameter which absorbs the free parameter AC 2 of Ref. [35] . 
We can alter the behavior of the coupling, and thus the energy dependence of the solution, by adjusting the value of 
Aq CD . According to equation (17), the behavior of the saturation scale asymptotes to 



where \ T 



I 2x(7o) 

0(1 -7b)' 



(20) 



This is qualitatively different from the large- Y" behavior of the GBW model, thus it prevents us from reproducing the 
asymptotic behavior of that model. Nevertheless, we can choose a value of Aqcd that approximately matches to 
its value from the GBW model over the upper range of rapidities we calculate, 10 < Y < 20. Some experimentation 
: 0.001 GeV 2 gives a reasonable match. The corresponding value from the fits in J35J[3H] would be 

( °' 24 4(6) CV) ' = 0.002 GeV 2 , which is fairly close. 
Note that the value of the effective Aqcd parameter determined by our fit is fairly small, which possibly indicates 
the need for further higher order corrections in the BK evolution (apart from the running coupling). 

The solution we have obtained for the BK evolution with running coupling is shown in figure [4] Although the 
overall shapes of the curves are qualitatively similar, the evolution with rapidity is slowed by the running coupling. 
Overall, the actual value of F x is less at large values of k 2 than it was in the fixed coupling case. 



shows that Aq CD 
A 2 
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FIG. 4: The gluon distribution from the numerical BK evolution with the running coupling for central p-Pb collisions. This is 
a continuation of figure [2j and the labels have the same interpretation. As is well known, the running coupling slows down the 
evolution of the wavefront. Otherwise, the plots are similar to the fixed-coupling BK evolution, though there is some erratic 
behavior at lower momenta especially at high rapidities which was not seen to such a large extent in the fixed coupling case. 



IV. RESULTS 



With the gluon distribution determined from either the analytic GBW model and the numerical solution of the BK 



evolution described in section III we can compute and study the correlation defined in section [HI There are eight 



parameters involved in this calculation: 

• the virtual photon mass M 

• the photon and pion rapidities Y 1 and Y n 

• the profile function c{b) 

• the mass number A 

• the center-of-mass energy y^i per nucleon 

• the detector-imposed cutoff on the transverse momenta of the virtual photons or pions, pj_ C ut 

• Different types of projectiles: either deuterons or protons. 

To simulate results from d-Au collisions at RHIC, we set s/s = 200 GeV, A = 197, and p± cu t — 1.5 GeV together 
with the use of the deuteron parton distribution functions. For LHC p-Pb collisions, we put y/s — 8800 GeV which 
corresponds to the full design energy (14TeV pp CM energy), A = 208, and pj_ cut = 3 GeV. For both the RHIC 
and LHC measurement, we compute the correlation for a low and a high value of M. In addition, we set Y" 7 = Y n 
which means that the virtual photons and pions are produced in the same rapidity window, and use c(b) = 0.85 which 
corresponds to central collisions. 

We show our final results for the correlations in figures [5] and [6] with different specifically chosen parameters. From 
figs. [5] and [6j we find that indeed the GBW model and the BK equation with fixed and running couplings yield similar 
double-peak structure, as expected from our early discussion for the correlation function at the away side (Acj) — n). 
However, they differ dramatically at the near side (A<p ~ or 27r). It is easy to see that the near side correlation 
samples the large q± region of the gluon distribution F Xg (q±). Since this distribution is slightly greater with the 
running coupling than with fixed coupling, the near side peak is somewhat enhanced in the running coupling case. 

We believe that the solution to the BK equation gives the correct large q± behavior for the dipole gluon distribution 
while the GBW model fails at large q±. Therefore, we expect that the curves generated by the numeric solution of 
the BK equation predicts the correlation for RHIC and LHC measurement for all azimuthal angle including both the 
near side and the away side. 
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V. CONCLUSION 

In this paper, by using the GBW model and the numerical solution to the BK equation, we have numerically studied 
the correlations between the lepton-pair and associated hadron in Drcll-Yan processes in pA collisions, which directly 
probe the dipole gluon distributions. The correlations have been computed for both RHIC d-Au collisions and LHC 
p-Pb collisions. We find that the Drell-Yan type correlation implies a strong suppression of the away side hadron 
at forward rapidity due to strong interaction between the forward quark from the projectile proton and the dense 
gluon from the target nucleus. In addition, we also emphasize that the Drell-Yan type correlation exhibits a unique 
double-peak structure in the away side. 
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FIG. 5: The angular correlations between the virtual photons and pions at RHIC, at medium rapidity, Yy = Y n = 2.5. The 
upper graph shows the correlation for a virtual photon mass of M = 0.5 GeV, and the lower one, for M = 4GeV. In each 
case, the three curves for GBW, fixed coupling BK, and running coupling BK, exhibit basically the same double-peak structure 
around A(f> = it, but they show differing behavior near A(j> = 0, 2ir, the near side correlation. This relates to the large-fc 2 
behavior of the corresponding gluon distributions. 




FIG. 6: The virtual photon-pion angular correlations at LHC at rapidity Yy = Y v = 4. The upper and lower graphs show 
M = 4 GeV and M = 8 GeV, respectively. As in figure [5] the key difference is in the near side correlation, arising due to 
differences in the high-momentum region of the gluon distribution. 



